# Essential libraries
library(corrplot)
library(ggeffects)
library(ggplot2)
library(DHARMa)
library(glmmTMB)
library(dplyr)
library(lubridate)
library(emmeans)
library(tidyverse)
library(performance)
df <- read.csv("./data/Mock_Mosquito_MixedEffects_w_slopes_Dataset.csv")
df$Trap_type = factor(df$Trap_type)
df$Mosquito_species = factor(df$Mosquito_species)
head(df)
We start by first defining the random effects part. There are two random effects defined in our model, one for the nested structure in our data and the other for the mosquito species. By defining random intercepts in m1 model, we assume that different cities have different baseline abundance levels, different sites within the same city have different baseline abundance, and different mosquito species have different baseline abundance levels. However, the slopes remain consistent meaning all cities and species respond to green coverage the same way.
m1_baseline <- glmmTMB(
Mosquito_abundance ~ Green_area_coverage * Trap_type +
Temperature * Relative_humidity +
Host_abundance * Predator_abundance +
Blue_area_coverage +
factor(Month) +
(1 | City / Location.Site) + # nested random effects
(1 | Mosquito_species), # random intercept by species
family = nbinom2(),
data = df,
REML = TRUE
)
summary(m1_baseline)
Family: nbinom2 ( log )
Formula: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature *
Relative_humidity + Host_abundance * Predator_abundance +
Blue_area_coverage + factor(Month) + (1 | City/Location.Site) + (1 | Mosquito_species)
Data: df
AIC BIC logLik -2*log(L) df.resid
18366.2 18515.3 -9157.1 18314.2 2262
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
Location.Site:City (Intercept) 0.12116 0.3481
City (Intercept) 0.12336 0.3512
Mosquito_species (Intercept) 0.05108 0.2260
Number of obs: 2288, groups: Location.Site:City, 54; City, 6; Mosquito_species, 5
Dispersion parameter for nbinom2 family (): 1.98
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.867e+00 7.967e-01 -2.344 0.019084 *
Green_area_coverage 8.413e-03 2.284e-03 3.683 0.000231 ***
Trap_typeLight_CO2 1.175e-01 7.294e-02 1.611 0.107204
Temperature 7.078e-02 3.399e-02 2.082 0.037300 *
Relative_humidity 3.435e-02 8.128e-03 4.226 2.38e-05 ***
Host_abundance 2.017e-02 9.056e-03 2.227 0.025963 *
Predator_abundance 2.250e-03 6.004e-03 0.375 0.707912
Blue_area_coverage 6.065e-03 3.833e-03 1.583 0.113513
factor(Month)2 -5.252e-02 8.426e-02 -0.623 0.533063
factor(Month)3 1.012e-01 9.255e-02 1.093 0.274366
factor(Month)4 1.588e-01 1.064e-01 1.492 0.135757
factor(Month)5 3.574e-01 1.201e-01 2.975 0.002934 **
factor(Month)6 2.764e-01 1.308e-01 2.113 0.034565 *
factor(Month)7 3.973e-01 1.215e-01 3.270 0.001075 **
factor(Month)8 1.085e-01 1.059e-01 1.025 0.305318
factor(Month)9 1.525e-01 9.086e-02 1.678 0.093279 .
factor(Month)10 8.593e-02 8.213e-02 1.046 0.295435
factor(Month)11 3.149e-04 7.928e-02 0.004 0.996830
factor(Month)12 -1.800e-03 8.171e-02 -0.022 0.982421
Green_area_coverage:Trap_typeLight_CO2 5.819e-03 1.265e-03 4.599 4.26e-06 ***
Temperature:Relative_humidity -8.019e-04 4.515e-04 -1.776 0.075746 .
Host_abundance:Predator_abundance 1.507e-05 1.160e-04 0.130 0.896639
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
baseline_aic <- AIC(m1_baseline)
cat("Baseline AIC:", baseline_aic, "\n")
Baseline AIC: 18366.19
In the second model, we consider random slopes for green levels in the cities, meaning that the effect of urban green differs per city.
m2_city_green <- glmmTMB(
Mosquito_abundance ~ Green_area_coverage * Trap_type +
Temperature * Relative_humidity +
Host_abundance * Predator_abundance +
Blue_area_coverage +
factor(Month) +
(1 + Green_area_coverage | City) + # City random intercept + slope
(1 | Location.Site:City) +
(1 | Mosquito_species),
family = nbinom2(),
data = df,
REML =TRUE
)
# Check convergence
check_convergence(m2_city_green)
[1] TRUE
m2_city_green$sdr$pdHess
[1] TRUE
cat("AIC difference:", AIC(m2_city_green) - baseline_aic, "\n")
AIC difference: -2.716824
lrt_m2 <- anova(m1_baseline, m2_city_green)$Chisq[2]
p_m2 <- 0.5 * pchisq(lrt_m2, df = 2, lower.tail = FALSE) +
0.5 * pchisq(lrt_m2, df = 1, lower.tail = FALSE)
cat("p-value:", p_m2, "\n")
p-value: 0.02217061
AIC difference is negative and the \(H_0: \text{Cities do not differ in how urban green level affects mosquito abundance}\) is rejected (p<0.05), the random slope of urban green level is significant.
m3_city_temp <- glmmTMB(
Mosquito_abundance ~ Green_area_coverage * Trap_type +
Temperature * Relative_humidity +
Host_abundance * Predator_abundance +
Blue_area_coverage +
factor(Month) +
(1 + Temperature | City) +
(1 | Location.Site:City) +
(1 | Mosquito_species),
family = nbinom2(),
data = df,
REML = TRUE
)
# Check convergence
check_convergence(m3_city_temp)
[1] FALSE
m3_city_temp$sdr$pdHess
[1] TRUE
cat("AIC difference:", AIC(m3_city_temp) - baseline_aic, "\n")
AIC difference: 3.928898
lrt_m3 <- anova(m1_baseline, m3_city_temp)$Chisq[2]
p_m3 <- 0.5 * pchisq(lrt_m3, df = 2, lower.tail = FALSE) +
0.5 * pchisq(lrt_m3, df = 1, lower.tail = FALSE)
cat("p-value for Temp slope:", p_m3, "\n")
p-value for Temp slope: 0.8774065
AIC difference is positive and the \(H_0: \text{Cities do not differ in how temperature affects mosquito abundance}\) is retained (p>0.05), the random slope of temperature is not significant.
m4_both_slopes <- glmmTMB(
Mosquito_abundance ~ Green_area_coverage * Trap_type +
Temperature * Relative_humidity +
Host_abundance * Predator_abundance +
Blue_area_coverage +
factor(Month) +
(1 + Green_area_coverage + Temperature | City) +
(1 | Location.Site:City) +
(1 | Mosquito_species),
family = nbinom2(),
data = df,
REML = TRUE
)
# Check convergence
check_convergence(m4_both_slopes)
[1] FALSE
m4_both_slopes$sdr$pdHess
[1] TRUE
# Compare m2 -> m4 (does adding temperature slope to m2 improve fit?)
lrt_m2_m4 <- anova(m2_city_green, m4_both_slopes)$Chisq[2]
cat("LRT stat (m2 -> m4):", lrt_m2_m4, "\n")
LRT stat (m2 -> m4): 0.204401
p_m2_m4 <- 0.5 * pchisq(lrt_m2_m4, df = 2, lower.tail = FALSE) +
0.5 * pchisq(lrt_m2_m4, df = 1, lower.tail = FALSE)
cat("p-value (m2 -> m4):", p_m2_m4, "\n")
p-value (m2 -> m4): 0.7770201
# Compare m3 -> m4 (does adding Green slope to m3 improve fit?)
lrt_m3_m4 <- anova(m3_city_temp, m4_both_slopes)$Chisq[2]
cat("LRT stat (m3 -> m4):", lrt_m3_m4, "\n")
LRT stat (m3 -> m4): 6.850123
p_m3_m4 <- 0.5 * pchisq(lrt_m3_m4, df = 2, lower.tail = FALSE) +
0.5 * pchisq(lrt_m3_m4, df = 1, lower.tail = FALSE)
cat("p-value (m3 -> m4):", p_m3_m4, "\n")
p-value (m3 -> m4): 0.0207054
# Compare all models by AIC
models <- list(m1_baseline = m1_baseline, m2 = m2_city_green, m3 = m3_city_temp, m4 = m4_both_slopes)
aic_tab <- data.frame(
model = names(models),
AIC = sapply(models, AIC),
pdHess = sapply(models, function(m) { if (is.null(m$sdr)) NA else m$sdr$pdHess })
)
aic_tab <- aic_tab %>% arrange(AIC)
print(aic_tab)
Across the models tested, the lowest AIC value was tested was achieved by model 2, for the remainder of the model selection procedure we will use the random effects defined in model 2. Now, we will focus on getting a more parsimonious fixed effects part by removing interaction terms in a step wise manner.
m_full <- glmmTMB(
Mosquito_abundance ~ Green_area_coverage * Trap_type +
Temperature * Relative_humidity +
Host_abundance * Predator_abundance +
Blue_area_coverage +
factor(Month) +
(1 + Green_area_coverage | City) + # City random intercept + slope
(1 | Location.Site:City) +
(1 | Mosquito_species),
family = nbinom2(),
data = df,
REML = FALSE # ML for fixed effects comparison
)
summary(m_full)
Family: nbinom2 ( log )
Formula: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature *
Relative_humidity + Host_abundance * Predator_abundance +
Blue_area_coverage + factor(Month) + (1 + Green_area_coverage |
City) + (1 | Location.Site:City) + (1 | Mosquito_species)
Data: df
AIC BIC logLik -2*log(L) df.resid
18216.1 18376.7 -9080.0 18160.1 2260
Random effects:
Conditional model:
Groups Name Variance Std.Dev. Corr
City (Intercept) 6.083e-02 0.246641
Green_area_coverage 3.707e-05 0.006088 -0.37
Location.Site:City (Intercept) 8.948e-02 0.299125
Mosquito_species (Intercept) 4.596e-02 0.214383
Number of obs: 2288, groups: City, 6; Location.Site:City, 54; Mosquito_species, 5
Dispersion parameter for nbinom2 family (): 2
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.862e+00 7.828e-01 -2.379 0.017359 *
Green_area_coverage 7.260e-03 3.286e-03 2.209 0.027153 *
Trap_typeLight_CO2 1.245e-01 7.270e-02 1.713 0.086802 .
Temperature 6.912e-02 3.385e-02 2.042 0.041182 *
Relative_humidity 3.430e-02 8.095e-03 4.237 2.27e-05 ***
Host_abundance 2.047e-02 8.994e-03 2.276 0.022868 *
Predator_abundance 2.477e-03 5.961e-03 0.415 0.677825
Blue_area_coverage 7.381e-03 3.584e-03 2.060 0.039431 *
factor(Month)2 -4.997e-02 8.398e-02 -0.595 0.551836
factor(Month)3 1.074e-01 9.227e-02 1.164 0.244285
factor(Month)4 1.671e-01 1.062e-01 1.574 0.115443
factor(Month)5 3.700e-01 1.198e-01 3.088 0.002016 **
factor(Month)6 2.886e-01 1.304e-01 2.213 0.026866 *
factor(Month)7 4.082e-01 1.212e-01 3.368 0.000757 ***
factor(Month)8 1.170e-01 1.056e-01 1.108 0.268040
factor(Month)9 1.577e-01 9.056e-02 1.741 0.081608 .
factor(Month)10 8.754e-02 8.185e-02 1.070 0.284832
factor(Month)11 5.986e-04 7.900e-02 0.008 0.993954
factor(Month)12 -1.046e-03 8.141e-02 -0.013 0.989752
Green_area_coverage:Trap_typeLight_CO2 5.691e-03 1.262e-03 4.511 6.45e-06 ***
Temperature:Relative_humidity -7.912e-04 4.496e-04 -1.760 0.078454 .
Host_abundance:Predator_abundance 1.126e-05 1.152e-04 0.098 0.922116
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
We start by removing the least significant interaction term which is Host_abundance:Predator_abundance .
m_wo_host_pre <- glmmTMB(
Mosquito_abundance ~ Green_area_coverage * Trap_type +
Temperature * Relative_humidity +
Host_abundance + Predator_abundance +
Blue_area_coverage +
factor(Month) +
(1 + Green_area_coverage | City) + # City random intercept + slope
(1 | Location.Site:City) +
(1 | Mosquito_species),
family = nbinom2(),
data = df,
REML = FALSE # ML for fixed effects comparison
)
summary(m_wo_host_pre)
Family: nbinom2 ( log )
Formula: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature *
Relative_humidity + Host_abundance + Predator_abundance +
Blue_area_coverage + factor(Month) + (1 + Green_area_coverage |
City) + (1 | Location.Site:City) + (1 | Mosquito_species)
Data: df
AIC BIC logLik -2*log(L) df.resid
18214.1 18369.0 -9080.0 18160.1 2261
Random effects:
Conditional model:
Groups Name Variance Std.Dev. Corr
City (Intercept) 6.066e-02 0.246290
Green_area_coverage 3.711e-05 0.006092 -0.37
Location.Site:City (Intercept) 8.961e-02 0.299345
Mosquito_species (Intercept) 4.597e-02 0.214411
Number of obs: 2288, groups: City, 6; Location.Site:City, 54; Mosquito_species, 5
Dispersion parameter for nbinom2 family (): 2
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.9033561 0.6709093 -2.837 0.004554 **
Green_area_coverage 0.0072419 0.0032831 2.206 0.027399 *
Trap_typeLight_CO2 0.1243781 0.0726983 1.711 0.087103 .
Temperature 0.0690642 0.0349110 1.978 0.047896 *
Relative_humidity 0.0342986 0.0082954 4.135 3.56e-05 ***
Host_abundance 0.0213223 0.0021267 10.026 < 2e-16 ***
Predator_abundance 0.0030338 0.0017445 1.739 0.082024 .
Blue_area_coverage 0.0074033 0.0035795 2.068 0.038618 *
factor(Month)2 -0.0501047 0.0839359 -0.597 0.550547
factor(Month)3 0.1074290 0.0921659 1.166 0.243775
factor(Month)4 0.1675267 0.1060089 1.580 0.114037
factor(Month)5 0.3701742 0.1199362 3.086 0.002026 **
factor(Month)6 0.2887538 0.1304774 2.213 0.026894 *
factor(Month)7 0.4082087 0.1213271 3.365 0.000767 ***
factor(Month)8 0.1170888 0.1055235 1.110 0.267172
factor(Month)9 0.1578000 0.0904607 1.744 0.081089 .
factor(Month)10 0.0876332 0.0818319 1.071 0.284218
factor(Month)11 0.0006265 0.0790036 0.008 0.993673
factor(Month)12 -0.0008575 0.0813868 -0.011 0.991593
Green_area_coverage:Trap_typeLight_CO2 0.0056932 0.0012615 4.513 6.38e-06 ***
Temperature:Relative_humidity -0.0007906 0.0004658 -1.697 0.089621 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
anova(m_full, m_wo_host_pre)
Data: df
Models:
m_wo_host_pre: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature * , zi=~0, disp=~1
m_wo_host_pre: Relative_humidity + Host_abundance + Predator_abundance + , zi=~0, disp=~1
m_wo_host_pre: Blue_area_coverage + factor(Month) + (1 + Green_area_coverage | , zi=~0, disp=~1
m_wo_host_pre: City) + (1 | Location.Site:City) + (1 | Mosquito_species), zi=~0, disp=~1
m_full: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature * , zi=~0, disp=~1
m_full: Relative_humidity + Host_abundance * Predator_abundance + , zi=~0, disp=~1
m_full: Blue_area_coverage + factor(Month) + (1 + Green_area_coverage | , zi=~0, disp=~1
m_full: City) + (1 | Location.Site:City) + (1 | Mosquito_species), zi=~0, disp=~1
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
m_wo_host_pre 27 18214 18369 -9080 18160
m_full 28 18216 18377 -9080 18160 0.0095 1 0.9222
Removing the Host_abundance:Predator_abundance interaction does not significantly worsen model fit, then the interaction is not needed. We move on to the next interaction.
m_wo_temp_hum <- glmmTMB(
Mosquito_abundance ~ Green_area_coverage * Trap_type +
Temperature + Relative_humidity +
Host_abundance + Predator_abundance +
Blue_area_coverage +
factor(Month) +
(1 + Green_area_coverage | City) + # City random intercept + slope
(1 | Location.Site:City) +
(1 | Mosquito_species),
family = nbinom2(),
data = df,
REML = FALSE # ML for fixed effects comparison
)
summary(m_wo_temp_hum)
Family: nbinom2 ( log )
Formula: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature +
Relative_humidity + Host_abundance + Predator_abundance +
Blue_area_coverage + factor(Month) + (1 + Green_area_coverage |
City) + (1 | Location.Site:City) + (1 | Mosquito_species)
Data: df
AIC BIC logLik -2*log(L) df.resid
18215.2 18364.3 -9081.6 18163.2 2262
Random effects:
Conditional model:
Groups Name Variance Std.Dev. Corr
City (Intercept) 6.149e-02 0.247970
Green_area_coverage 3.759e-05 0.006131 -0.38
Location.Site:City (Intercept) 8.909e-02 0.298473
Mosquito_species (Intercept) 4.585e-02 0.214137
Number of obs: 2288, groups: City, 6; Location.Site:City, 54; Mosquito_species, 5
Dispersion parameter for nbinom2 family (): 2
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.9280810 0.3606140 -2.574 0.010064 *
Green_area_coverage 0.0071653 0.0032923 2.176 0.029527 *
Trap_typeLight_CO2 0.1242485 0.0727127 1.709 0.087496 .
Temperature 0.0113019 0.0082549 1.369 0.170962
Relative_humidity 0.0211927 0.0031634 6.699 2.09e-11 ***
Host_abundance 0.0212425 0.0021284 9.981 < 2e-16 ***
Predator_abundance 0.0031070 0.0017451 1.780 0.075007 .
Blue_area_coverage 0.0074403 0.0035686 2.085 0.037073 *
factor(Month)2 -0.0599509 0.0838062 -0.715 0.474392
factor(Month)3 0.0958339 0.0921129 1.040 0.298156
factor(Month)4 0.1601716 0.1060906 1.510 0.131104
factor(Month)5 0.3753467 0.1198993 3.131 0.001745 **
factor(Month)6 0.2967466 0.1303588 2.276 0.022823 *
factor(Month)7 0.4141824 0.1212072 3.417 0.000633 ***
factor(Month)8 0.1090555 0.1055911 1.033 0.301693
factor(Month)9 0.1453758 0.0903559 1.609 0.107633
factor(Month)10 0.0820182 0.0818185 1.002 0.316131
factor(Month)11 0.0004793 0.0790249 0.006 0.995160
factor(Month)12 -0.0005138 0.0814086 -0.006 0.994964
Green_area_coverage:Trap_typeLight_CO2 0.0057447 0.0012610 4.556 5.22e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
anova(m_wo_host_pre, m_wo_temp_hum)
Data: df
Models:
m_wo_temp_hum: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature + , zi=~0, disp=~1
m_wo_temp_hum: Relative_humidity + Host_abundance + Predator_abundance + , zi=~0, disp=~1
m_wo_temp_hum: Blue_area_coverage + factor(Month) + (1 + Green_area_coverage | , zi=~0, disp=~1
m_wo_temp_hum: City) + (1 | Location.Site:City) + (1 | Mosquito_species), zi=~0, disp=~1
m_wo_host_pre: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature * , zi=~0, disp=~1
m_wo_host_pre: Relative_humidity + Host_abundance + Predator_abundance + , zi=~0, disp=~1
m_wo_host_pre: Blue_area_coverage + factor(Month) + (1 + Green_area_coverage | , zi=~0, disp=~1
m_wo_host_pre: City) + (1 | Location.Site:City) + (1 | Mosquito_species), zi=~0, disp=~1
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
m_wo_temp_hum 26 18215 18364 -9081.6 18163
m_wo_host_pre 27 18214 18369 -9080.0 18160 3.1134 1 0.07765 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Removing the Temp x Humidity interaction does not worsen model fit, the interaction should be removed. The remaining interaction is significant, therefore we move to the main effects.
m_wo_blue <- glmmTMB(
Mosquito_abundance ~ Green_area_coverage * Trap_type +
Temperature + Relative_humidity +
Host_abundance + Predator_abundance +
factor(Month) +
(1 + Green_area_coverage | City) + # City random intercept + slope
(1 | Location.Site:City) +
(1 | Mosquito_species),
family = nbinom2(),
data = df,
REML = FALSE
)
summary(m_wo_blue)
Family: nbinom2 ( log )
Formula: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature +
Relative_humidity + Host_abundance + Predator_abundance +
factor(Month) + (1 + Green_area_coverage | City) + (1 | Location.Site:City) + (1 | Mosquito_species)
Data: df
AIC BIC logLik -2*log(L) df.resid
18217.4 18360.8 -9083.7 18167.4 2263
Random effects:
Conditional model:
Groups Name Variance Std.Dev. Corr
City (Intercept) 9.077e-02 0.30127
Green_area_coverage 3.587e-05 0.00599 -0.43
Location.Site:City (Intercept) 9.713e-02 0.31166
Mosquito_species (Intercept) 4.650e-02 0.21565
Number of obs: 2288, groups: City, 6; Location.Site:City, 54; Mosquito_species, 5
Dispersion parameter for nbinom2 family (): 2
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.7160550 0.3541022 -2.022 0.04316 *
Green_area_coverage 0.0074203 0.0032834 2.260 0.02382 *
Trap_typeLight_CO2 0.1237693 0.0726909 1.703 0.08863 .
Temperature 0.0112996 0.0082553 1.369 0.17107
Relative_humidity 0.0212518 0.0031635 6.718 1.85e-11 ***
Host_abundance 0.0211876 0.0021288 9.953 < 2e-16 ***
Predator_abundance 0.0030344 0.0017452 1.739 0.08208 .
factor(Month)2 -0.0604923 0.0838047 -0.722 0.47040
factor(Month)3 0.0956423 0.0921279 1.038 0.29920
factor(Month)4 0.1595994 0.1060774 1.505 0.13244
factor(Month)5 0.3750395 0.1198909 3.128 0.00176 **
factor(Month)6 0.2964331 0.1303665 2.274 0.02298 *
factor(Month)7 0.4142905 0.1211931 3.418 0.00063 ***
factor(Month)8 0.1085706 0.1055968 1.028 0.30387
factor(Month)9 0.1443938 0.0903574 1.598 0.11004
factor(Month)10 0.0821350 0.0818163 1.004 0.31543
factor(Month)11 0.0006321 0.0790305 0.008 0.99362
factor(Month)12 -0.0016955 0.0814145 -0.021 0.98339
Green_area_coverage:Trap_typeLight_CO2 0.0057879 0.0012602 4.593 4.37e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
anova(m_wo_temp_hum, m_wo_blue)
Data: df
Models:
m_wo_blue: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature + , zi=~0, disp=~1
m_wo_blue: Relative_humidity + Host_abundance + Predator_abundance + , zi=~0, disp=~1
m_wo_blue: factor(Month) + (1 + Green_area_coverage | City) + (1 | Location.Site:City) + , zi=~0, disp=~1
m_wo_blue: (1 | Mosquito_species), zi=~0, disp=~1
m_wo_temp_hum: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature + , zi=~0, disp=~1
m_wo_temp_hum: Relative_humidity + Host_abundance + Predator_abundance + , zi=~0, disp=~1
m_wo_temp_hum: Blue_area_coverage + factor(Month) + (1 + Green_area_coverage | , zi=~0, disp=~1
m_wo_temp_hum: City) + (1 | Location.Site:City) + (1 | Mosquito_species), zi=~0, disp=~1
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
m_wo_blue 25 18217 18361 -9083.7 18167
m_wo_temp_hum 26 18215 18364 -9081.6 18163 4.2241 1 0.03985 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Removal of blue zone coverage worsened the model, we retain it. Month shows mixed significance, so we will try removing it and compare the model fits.
m_wo_month <- glmmTMB(
Mosquito_abundance ~ Green_area_coverage * Trap_type +
Temperature + Relative_humidity +
Host_abundance + Predator_abundance +
Blue_area_coverage +
(1 + Green_area_coverage | City) + # City random intercept + slope
(1 | Location.Site:City) +
(1 | Mosquito_species),
family = nbinom2(),
data = df,
REML = FALSE
)
anova(m_wo_temp_hum, m_wo_month)
Data: df
Models:
m_wo_month: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature + , zi=~0, disp=~1
m_wo_month: Relative_humidity + Host_abundance + Predator_abundance + , zi=~0, disp=~1
m_wo_month: Blue_area_coverage + (1 + Green_area_coverage | City) + (1 | , zi=~0, disp=~1
m_wo_month: Location.Site:City) + (1 | Mosquito_species), zi=~0, disp=~1
m_wo_temp_hum: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature + , zi=~0, disp=~1
m_wo_temp_hum: Relative_humidity + Host_abundance + Predator_abundance + , zi=~0, disp=~1
m_wo_temp_hum: Blue_area_coverage + factor(Month) + (1 + Green_area_coverage | , zi=~0, disp=~1
m_wo_temp_hum: City) + (1 | Location.Site:City) + (1 | Mosquito_species), zi=~0, disp=~1
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
m_wo_month 15 18223 18309 -9096.4 18193
m_wo_temp_hum 26 18215 18364 -9081.6 18163 29.623 11 0.001815 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Removing the month worsens the fit, so our final model is:
m_final <- glmmTMB(
Mosquito_abundance ~ Green_area_coverage * Trap_type +
Temperature + Relative_humidity +
Host_abundance + Predator_abundance +
factor(Month) +
Blue_area_coverage +
(1 + Green_area_coverage | City) + # City random intercept + slope
(1 | Location.Site:City) +
(1 | Mosquito_species),
family = nbinom2(),
data = df,
REML = TRUE
)
summary(m_final)
Family: nbinom2 ( log )
Formula: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature +
Relative_humidity + Host_abundance + Predator_abundance +
factor(Month) + Blue_area_coverage + (1 + Green_area_coverage |
City) + (1 | Location.Site:City) + (1 | Mosquito_species)
Data: df
AIC BIC logLik -2*log(L) df.resid
18332.7 18481.8 -9140.4 18280.7 2262
Random effects:
Conditional model:
Groups Name Variance Std.Dev. Corr
City (Intercept) 8.433e-02 0.290393
Green_area_coverage 5.007e-05 0.007076 -0.47
Location.Site:City (Intercept) 9.124e-02 0.302060
Mosquito_species (Intercept) 5.091e-02 0.225636
Number of obs: 2288, groups: City, 6; Location.Site:City, 54; Mosquito_species, 5
Dispersion parameter for nbinom2 family (): 1.98
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.9306526 0.3687708 -2.524 0.011614 *
Green_area_coverage 0.0070819 0.0036047 1.965 0.049460 *
Trap_typeLight_CO2 0.1241486 0.0729250 1.702 0.088677 .
Temperature 0.0116819 0.0082860 1.410 0.158586
Relative_humidity 0.0212013 0.0031738 6.680 2.39e-11 ***
Host_abundance 0.0212067 0.0021355 9.930 < 2e-16 ***
Predator_abundance 0.0030975 0.0017508 1.769 0.076866 .
factor(Month)2 -0.0609114 0.0840613 -0.725 0.468693
factor(Month)3 0.0934745 0.0924016 1.012 0.311724
factor(Month)4 0.1568507 0.1064307 1.474 0.140553
factor(Month)5 0.3708248 0.1202928 3.083 0.002051 **
factor(Month)6 0.2913307 0.1308038 2.227 0.025932 *
factor(Month)7 0.4096605 0.1216106 3.369 0.000755 ***
factor(Month)8 0.1055351 0.1059346 0.996 0.319139
factor(Month)9 0.1430460 0.0906455 1.578 0.114547
factor(Month)10 0.0811003 0.0820680 0.988 0.323051
factor(Month)11 0.0003518 0.0792615 0.004 0.996458
factor(Month)12 -0.0002625 0.0816592 -0.003 0.997435
Blue_area_coverage 0.0073499 0.0036179 2.032 0.042202 *
Green_area_coverage:Trap_typeLight_CO2 0.0057406 0.0012649 4.538 5.67e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
library(DHARMa)
# Simulate residuals
sim_res <- simulateResiduals(m_final)
# Plot residual diagnostics
plot(sim_res)
testDispersion(sim_res)
DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
data: simulationOutput
dispersion = 0.67548, p-value = 0.728
alternative hypothesis: two.sided
testZeroInflation(sim_res)
DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted
model
data: simulationOutput
ratioObsSim = 1.0195, p-value = 0.912
alternative hypothesis: two.sided
ranef(m_final)
$City
(Intercept) Green_area_coverage
Barcelona 0.05665254 -0.0006839102
Casablanca -0.23987633 -0.0074754880
Montpellier -0.28248408 0.0087906133
Nijmegen 0.30877301 -0.0061932812
Rabat 0.06450358 0.0013067294
Utrecht 0.09243128 0.0042553368
$`Location.Site:City`
(Intercept)
1:Barcelona 0.07043914
1:Casablanca 0.03530279
1:Montpellier -0.14612773
1:Nijmegen 0.01928858
1:Rabat 0.17883180
1:Utrecht -0.14043348
2:Barcelona -0.21674606
2:Casablanca 0.22542755
2:Montpellier -0.34050646
2:Nijmegen 0.52709126
2:Rabat 0.47659812
2:Utrecht 0.01212399
3:Barcelona 0.23358215
3:Casablanca -0.17418438
3:Montpellier -0.01351448
3:Nijmegen -0.58369434
3:Rabat -0.06142319
3:Utrecht 0.17008047
4:Barcelona -0.04615371
4:Casablanca -0.01484421
4:Montpellier 0.30291202
4:Nijmegen 0.34686847
4:Rabat -0.62301501
4:Utrecht -0.05344528
5:Barcelona 0.09420243
5:Casablanca -0.42233337
5:Montpellier -0.06700641
5:Nijmegen 0.18584290
5:Rabat 0.15272034
5:Utrecht 0.20916035
6:Barcelona -0.05731574
6:Casablanca 0.16563333
6:Montpellier 0.23113679
6:Nijmegen -0.27094925
6:Rabat 0.02313622
6:Utrecht -0.13363323
7:Barcelona -0.18372403
7:Casablanca -0.48904874
7:Montpellier -0.34834010
7:Nijmegen -0.02084802
7:Rabat -0.31336564
7:Utrecht -0.01039561
8:Barcelona 0.20354565
8:Casablanca 0.21546823
8:Montpellier -0.15625847
8:Nijmegen 0.22299093
8:Rabat 0.47731495
8:Utrecht 0.13005884
9:Barcelona -0.03759367
9:Casablanca -0.07094313
9:Montpellier 0.37974388
9:Nijmegen -0.16354145
9:Rabat -0.18702398
9:Utrecht 0.05690802
$Mosquito_species
(Intercept)
Aedes aegypti 0.38422767
Aedes albopictus -0.14656346
Anopheles gambiae -0.04100560
Culex pipiens -0.16165110
Culex quinquefasciatus -0.03500751
library(ggeffects)
preds <- ggpredict(m_final, terms = c("Green_area_coverage", "Trap_type"))
plot(preds) + labs(y="Predicted Mosquito Abundance")
pred_temp <- ggpredict(m_final, terms = "Temperature [all]")
plot(pred_temp) +
labs(title = "Predicted Mosquito Abundance",
y="Predicted Mosquito Abundance")
pred_host <- ggpredict(m_final, terms = "Host_abundance")
plot(pred_host) +
labs(x="Host abundance", y="Predicted mosquitoes") +
theme_minimal()
pred_pre <- ggpredict(m_final, terms = "Predator_abundance")
plot(pred_pre) +
labs(x="Predator abundance", y="Predicted mosquitoes") +
theme_minimal()
pred_month <- ggpredict(m_final, terms = "Month [1:12]")
plot(pred_month) +
scale_x_continuous(breaks = 1:12,
labels = month.abb) +
labs(x="Month", y="Predicted mosquitoes",
title="Seasonal effect") +
theme_minimal()
NA
NA
re_plots = sjPlot::plot_model(m_final, type = "re")
re_plots[[1]] + theme_minimal()
re_plots[[3]] + theme_minimal()
library(sjPlot)
plot_model(m_final, type = "est",
axis.lim = c(0.8, 1.4),
show.values = TRUE,
value.offset = 0.3,
value.size = 2.5,
dot.size = 2,
line.size = 0.7,
vline.color = 'black') +
labs(title = "Incidence Rate Ratios (IRR)",
x = "Incidence Rate Ratio") +
theme_minimal()
icc_values <- performance::icc(m_final, tolerance = 1e-6)
icc_values
# Intraclass Correlation Coefficient
Adjusted ICC: 0.399
Unadjusted ICC: 0.261
ggplot(df, aes(x = factor(Green_level), y = Mosquito_abundance, fill = Trap_type)) +
geom_boxplot(alpha = 0.7, outlier.size = 0.5) +
theme_minimal(base_size = 14) +
labs(
title = "Observed Mosquito Count by Greenness Level and Trap Type",
x = "Urban Greenness Level",
y = "Observed Mosquito Abundance"
) +
scale_fill_brewer(palette = "Set2")
Mosquito abundance increases consistently with greenness level across both trap types. Although Light–CO₂ traps and Human–Odor–CO₂ traps show similar distributions, greener areas clearly host higher mosquito numbers. This provides initial, descriptive evidence of an ecological gradient before modeling.
num_df <- df %>% select(Temperature_2m, Relative_humidity, Blue_area_coverage, Green_area_coverage = Green_level)
corrplot(cor(num_df), method = "color", type = "upper", tl.col = "black")
df$Green_level <- as.factor(df$Green_level)
# Refit model
model_abundance2 <- glmmTMB(
Mosquito_abundance ~ Green_level * Trap_type + Temperature_2m +
Relative_humidity + Blue_area_coverage +
(1 | Location.Site) + (1 | Mosquito_species),
family = nbinom2(),
data = df
)
# Estimated marginal means
emm <- emmeans(model_abundance2, ~ Green_level * Trap_type, type = "response")
# Plot interaction
plot(emm, comparison = TRUE) +
labs(
title = "Predicted Mosquito Abundance by Green Level and Trap Type",
x = "Urban Greenness Level",
y = "Predicted Abundance (Response Scale)"
)
The model predicts a strong increase in mosquito abundance with higher urban greenness, indicating vegetation-rich environments favor mosquito populations. The effect of trap type is minimal — both traps capture similar numbers of mosquitoes across all vegetation levels. Confidence intervals widen at higher greenness levels, suggesting greater variability in mosquito counts in green zones.
# Extract random effects and standard deviations
ran_list <- ranef(model_abundance)$cond
# Convert to a tidy data frame
ran_df <- bind_rows(
lapply(names(ran_list), function(g) {
data.frame(
Group = g,
Level = rownames(ran_list[[g]]),
Estimate = ran_list[[g]][, "(Intercept)"]
)
})
)
# Plot
ggplot(ran_df, aes(x = reorder(Level, Estimate), y = Estimate, color = Group)) +
geom_point(size = 2) +
coord_flip() +
facet_wrap(~ Group, scales = "free_y") +
theme_minimal(base_size = 13) +
labs(
title = "Random Effects by Site and Species",
x = "Level (Site / Species)",
y = "Deviation from Average Abundance"
) +
theme(legend.position = "none",
axis.text.y = element_text(angle = 15, vjust = 0.5),)
ggplot(df, aes(x = Temperature_2m)) +
geom_histogram(binwidth = 1, fill = "steelblue", color = "white", alpha = 0.8) +
labs(
title = "Temperature Distribution Across Sampling Sites",
x = "Temperature (°C)",
y = "Frequency"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", size = 16),
axis.text.x = element_text(angle = 45, hjust = 1)
)
ggplot(df, aes(x = City, y = Relative_humidity, fill = City)) +
geom_boxplot(alpha = 0.7) +
labs(
title = "Variation in Relative Humidity by City",
x = "City",
y = "Relative Humidity (%)"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", size = 16),
legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1)
)
df %>%
ggplot(aes(x = Wind_speed_10m, fill = City)) +
geom_density(alpha = 0.5) +
labs(
title = "Distribution of Wind Speed Across Cities",
x = "Wind Speed (m/s)",
y = "Density",
fill = "City"
) +
theme_minimal(base_size = 14)
df %>%
mutate(Month = month(Sampling_date, label = TRUE, abbr = TRUE)) %>%
group_by(City, Month) %>%
summarise(Mean_Temperature = mean(Temperature_2m, na.rm = TRUE)) %>%
ggplot(aes(x = Month, y = Mean_Temperature, color = City, group = City)) +
geom_line(size = 1.2) +
geom_point(size = 2) +
labs(
title = "Average Monthly Temperature by City",
x = "Month",
y = "Temperature (°C)",
color = "City"
) +
theme_minimal(base_size = 14) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(face = "bold", size = 16)
)
pred_inter <- ggpredict(
model_abundance,
terms = c("Green_level", "Trap_type"),
type = "fixed", # correct for your version
bias_correction = TRUE
)
pred_df <- as.data.frame(pred_inter)
ggplot(pred_df, aes(
x = factor(x), # Greenness levels
y = group, # Trap types
fill = predicted
)) +
geom_tile(color = "white") +
geom_text(aes(label = round(predicted, 1)), color = "black", size = 4) +
scale_fill_gradient(low = "#fefefe", high = "#01665e") +
labs(
title = "Predicted Mosquito Count by Greenness Level and Trap Type",
x = "Urban Greenness Level",
y = "Trap Type",
fill = "Predicted\nAbundance"
) +
theme_minimal(base_size = 14)+
theme(plot.title = element_text(size=12))
Mosquito abundance rises markedly with increasing urban greenness, regardless of the trap type used. The model predicts roughly five times more mosquitoes in highly vegetated areas than in low-greenness zones, highlighting the ecological importance of vegetation in mosquito distribution.
df_pred <- df %>%
mutate(
Predicted = predict(model_abundance, type = "response"),
Observed = Mosquito_abundance
)
ggplot(df_pred, aes(x = Predicted, y = Observed)) +
geom_point(alpha = 0.4, size = 2) +
geom_smooth(method = "lm", color = "#01665e", se = TRUE) +
labs(
title = "Model Fit: Predicted vs Observed Mosquito Abundance",
x = "Predicted Abundance (Model)",
y = "Observed Abundance (Data)"
) +
theme_minimal(base_size = 14)
df %>%
group_by(Green_level, Mosquito_species) %>%
summarise(Mean_abundance = mean(Mosquito_abundance)) %>%
ggplot(aes(x = factor(Green_level), y = Mean_abundance, fill = Mosquito_species)) +
geom_bar(stat = "identity", position = "fill") +
labs(
title = "Relative Mosquito Species Composition by Greenness Level",
x = "Urban Greenness Level",
y = "Proportion of Total Abundance"
) +
theme_minimal(base_size = 14)
df %>%
mutate(Month = month(Sampling_date, label = TRUE, abbr = TRUE)) %>%
group_by(Month) %>%
summarise(Mean_abundance = mean(Mosquito_abundance)) %>%
ggplot(aes(x = Month, y = Mean_abundance, group = 1)) +
geom_line(size = 1.2, color = "#01665e") +
geom_point(size = 2, color = "#01665e") +
labs(
title = "Seasonal Trend in Mosquito Abundance",
x = "Month",
y = "Average Abundance"
) +
theme_minimal(base_size = 14)
df %>%
mutate(Month = month(Sampling_date, label = TRUE, abbr = TRUE)) %>%
group_by(City, Month) %>%
summarise(Mean_abundance = mean(Mosquito_abundance, na.rm = TRUE)) %>%
ggplot(aes(x = Month, y = Mean_abundance, group = 1)) +
geom_line(color = "#01665e", size = 1.2) +
geom_point(color = "#01665e", size = 2) +
facet_wrap(~ City, ncol = 2) +
labs(
title = "Seasonal Trends in Mosquito Abundance by City",
x = "Month",
y = "Average Mosquito Abundance"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", size = 16),
strip.text = element_text(size = 12),
axis.text.x = element_text(angle = 90, vjust = 0.5)
)
df %>%
mutate(Month = month(Sampling_date, label = TRUE, abbr = TRUE)) %>%
group_by(City, Month) %>%
summarise(Mean_abundance = mean(Mosquito_abundance, na.rm = TRUE)) %>%
ggplot(aes(x = Month, y = Mean_abundance, color = City, group = City)) +
geom_line(size = 1.2) +
geom_point(size = 2) +
labs(
title = "Seasonal Trends in Mosquito Abundance by City",
x = "Month",
y = "Average Mosquito Abundance",
color = "City"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", size = 16),
legend.position = "bottom"
)
df %>%
mutate(Month = month(Sampling_date, label = TRUE, abbr = TRUE)) %>%
group_by(Mosquito_species, Month) %>%
summarise(Mean_abundance = mean(Mosquito_abundance, na.rm = TRUE), .groups = "drop") %>%
ggplot(aes(x = Month, y = Mean_abundance, color = Mosquito_species, group = Mosquito_species)) +
geom_point(size = 2) +
geom_smooth(se = FALSE, linewidth = 1.2, method = "loess") +
labs(
title = "Seasonal Trends in Mosquito Abundance by Species",
x = "Month",
y = "Average Mosquito Abundance",
color = "Species"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", size = 16),
legend.position = "bottom"
)
df$Date <- as.Date(df$Sampling_date, format = "%Y-%m-%d")
# Create the monthly trend plot
df %>%
mutate(Month = month(Date, label = TRUE, abbr = TRUE)) %>%
group_by(Mosquito_species, Month) %>%
summarise(Mean_abundance = mean(Mosquito_abundance, na.rm = TRUE), .groups = "drop") %>%
ggplot(aes(x = Month, y = Mean_abundance, color = Mosquito_species, group = Mosquito_species)) +
geom_line(linewidth = 1.2) +
geom_point(size = 2) +
labs(
title = "Monthly Trends in Mosquito Abundance by Species",
x = "Month",
y = "Average Mosquito Abundance",
color = "Species"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", size = 16),
legend.position = "bottom"
)
# Create faceted plot by species
df %>%
mutate(Month = month(Date, label = TRUE, abbr = TRUE)) %>%
group_by(Mosquito_species, Month) %>%
summarise(Mean_abundance = mean(Mosquito_abundance, na.rm = TRUE), .groups = "drop") %>%
ggplot(aes(x = Month, y = Mean_abundance, group = 1)) +
geom_line(color = "#2C7BB6", linewidth = 1.2) +
geom_point(color = "#2C7BB6", size = 2) +
facet_wrap(~ Mosquito_species, scales = "free_y") +
labs(
title = "Monthly Trends in Mosquito Abundance by Species",
x = "Month",
y = "Average Mosquito Abundance"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
strip.text = element_text(face = "bold"),
axis.text.x = element_text(angle = 90, vjust = 0.5),
legend.position = "none"
)
ggplot(df, aes(x = Trap_type, y = Mosquito_abundance, fill = Trap_type)) +
geom_boxplot(outlier.alpha = 0.2) +
facet_wrap(~ City) +
scale_fill_brewer(palette = "Set2") +
labs(
title = "Distribution of Mosquito Abundance by Trap Type and City",
x = "Trap Type", y = "Mosquito Abundance"
) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1)
)